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ABSTRACT 

Radial and nonradial oscillations offer the opportunity to investigate the in- 
terior properties of stars. We use 2D stellar models and a 2D finite difference 
integration of the linearized pulsation equations to calculate non-radial oscilla- 
tions. This approach allows us to directly calculate the pulsation modes for a 
distorted rotating star without treating the rotation as a perturbation. We are 
also able to express the finite difference solution in the horizontal direction as a 
sum of multiple spherical harmonics for any given mode. Using these methods, 
we have investigated the effects of increasing rotation and the number of spherical 
harmonics on the calculated eigenfrequencies and eigenfunctions and compared 
the results to perturbation theory. In slowly rotating stars, current methods 
work well, and we show that the eigenfunction can be accurately modelled us- 
ing 2 nd order perturbation theory and a single spherical harmonic. We use 10 
M models with velocities ranging from to 420 km s _1 (0.89 Q c ) and examine 
low order p modes. We find that one spherical harmonic remains reasonable up 
to a rotation rate around 300km s _1 (0.69 Q c ) for the radial fundamental mode, 
but can fail at rotation rates as low as 90 km s _1 (0.23 Q c ) for the 2H mode or 
/ = 2 p 2 mode, based on the eigenfrequencies alone. Depending on the mode in 
question, a single spherical harmonic may fail at lower rotation rates if the shape 
of the eigenfunction is taken into consideration. Perturbation theory, in contrast, 
remains valid up to relatively high rotation rates for most modes. We find the 
lowest failure surface equatorial velocity is 120 km s _1 (0.30 Q c ) for the I = 2 p 2 
mode, but failure velocities between 240 and 300 km s _1 (0.58-0.69 f2 c )are more 
typical. 

Subject headings: stars: oscillations, stars: rotation 
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Introduction 



Stellar oscillations provide us with a probe of the internal structure of stars. The os- 
cillations depend on the stellar structure, and are modified by factors such as rotation, 
magnetic fields and tidal forces. In theory, if we have sufficiently accurate parameters for a 
star, we can produce models which will constrain the internal structure. Unfortunately, due 
to the uncertainties on the temperature and luminosity of the star and the large number 
of free parameters (mass, rotation rate, age, etc.), this process is much more difficult in 
practise. Accurate modelling also requires enough observed modes to actually place some 
constraints on the star. The more modes available, the tighter these constraints can be, 
but we must be sure that all the modes used are real. Artificial or extraneous modes can 
make it impossible to produce a matching model. In recent years, the number of stars with 
multiple modes has increased g reatly, both thanks to the grou nd based networks such as 



STEPHI ( IBelmonte et al. 



servations such as WIRE (Hacking et al. Ml999h and MOST flWalker et al. 1120031 ). Current 



1993 ]) and WET (jNather et al. II1990T). as well as space -based ob- 



and upcoming space m issions, such as Kepler (IBasri. Borucki fc Kochll2005l ) and COROT 



(IBaglin fc et al. 



20011 ) are expected to further increase the number of multiperiodic vari- 



ables. Unfortunately, the theory still lags behind the observations, particularly for rotating 
stars. 



The first investigation of non-radial oscillations was undertaken by IPekerid (119381 ). This 
paper derived the linearized, adiabatic equations for nonradial oscillations of non-rotating 
stars, and then solved the equations for models of uniform density. At the time, it was 
assumed that non-radial modes would be subject to significant amounts of damping, more 
so than the purely radial modes. As a result, non-radial oscillations were generally not 
studied extensively. However, these assumptions do not hold for the low order p modes or 
for all g modes. Unlike radial oscillations, which are unstable only for 7 < ~, there are 
some non-radial oscillations of a unifo r m den sity sphere which are unstable for all values 
of 7. Based on these results, IPekerid (119381) conclu ded that non-radial oscillations must 
be considered. Using these results, I Cowling (1l94ll ) calculated the periods of non-radial 
oscillations for non-rotating polytropes. 

Before the advent of numerical techniques, these equations had to be solved using ana- 
lytical methods. Much of this work was done by Chandrasekhar, who ex plored the variational 
princ i ple as a method of solv ing the linear adiabatic pulsation equations (I Chandrasekhar &: Lebovitz 
19621 ; I Chandrasekharl Il964l ) . This method depends on an arbitrary guess at the form of the 
eigenfunction, and the resulting eigenvalues depend on the guess. Fortunately, even marginal 
guesses at the eigenfunction can produce reasonable results for the eigenfrequencies with this 
method. This approach is largely unused today, as it has been superseded by computational 
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work using more efficient and accurate numerical techniques. 



The first direct num erical integration of the linearize d equations for nonradial oscil- 
lations was performed by iHurley. Roberts fc Wrightl (119661 ). In this work, they calculated 
oscillation frequencies for non-rotating, polytropic stellar models, for comparison with the 
earlier analytic approaches discussed above. Although they restricted themselves to poly- 
tropic models, their method can relatively easily be extended to more realistic stellar models. 

All of these approaches depend on perturbations to a non-rotating (i.e., spherical) stellar 
model. In this case, the calculations are relatively straight forward. Rotation, even moderate 
rotation, can significantly complicate the calculation, and many attempts have been made 
over the years to solve the problem with varying degrees of success. These will be discussed 
in more detail below. 

In spherical stars, the solution to the linear adiabatic pulsation equations is separable, 
and can be written as 

£ r = X(r)YT(M) (1) 

The angular variation can be characterized by a single spherical harmonic, Yf 1 , and both I 
and m are legitimate quantum numbers. Once a star becomes distorted, e.g. through tidal 
effects or rotation, the situation becomes more complex and several problems arise. The 
eigenfunction can no longer strictly be described by a single spherical harmonic, and thus / 
is no longer a valid quantum number. As long as the star remains axisymmetric, m remains 
valid. As well as changes in the structure of the eigenfunction, the pulsation frequencies 
themselves will change. It is this change in eigenfrequency that has been of most interest 
to researchers, particularly as observations continue to find more and more rotating and 
pulsating stars, many with multiple frequencies. 

One of the earli er attempts to solve the line a r adia batic pulsation equations for rotating 
stars was made by IChandrasekhar fc Lebovita (Il962l ). who applied the virial theorem to 
rotating incompressible fluids. T h e var i ationa l principle has also been extended to include 
slowly rotating stars by IClementl (11964 . If 9651 ) . Further attempt s at improving the method 
through a better choice of basis vectors have also been made by IClementl (119861 ) . Although 
the variational equations themselves can be applied to a star with an arbitrary rotation 
rate (fl), the method also depends on being able to model the structure of the star. The 
structure of rotating stars has generally been modelled as a perturbation to the non-rotating 
structure. Because the structural perturbations are limited to modelling slowly rotating 
stars, the variational method was also limited to slowly rotating stars. 

An approach us ed more frequently now is based on a perturbation approach, as devel- 
oped by ISaiol (ll98ll ). In this framework, the rotation is treated as a perturbation on the 
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structure of the star. For example, the radial location in a rotating model would be written 

as 

r = a[l + e(a, M)] (2) 

The linearized pulsation equations are expanded in a series in powers of the rotation rate. The 
zeroth power merely gives the nonrotating eigenvalues and eigenfunctions. Each non-rotating 
eigenfunction can be written in terms of a single spherical harmonic, and the eigenfunction 
can be characterized by three quantum numbers relating to the number of radial nodes and 
the two angular quantum numbers, / and m, associated with the spherical harmonics. The 
first order in the expansion in powers of the rotation rate lifts the 2/+1 fold degeneracy in 
the eigenvalues, while the eigenfunctions that correspond to this order are still characterized 
by a single spherical harmonic. 

We note that this will not be true in the general set of linearized pulsation equations of 
a rotating star. The coefficients of the perturbations in the pulsation equations, composed of 
terms based on the static rotating model, will have latitudinal variations. The eigenfunctions 
will also have a latitudinal variation, so that the equations can be expressed as products of 
spherical harmonics, which in turn can be written as sums of spherical harmonics through 
appropriate recursion relations. 

In perturbation theory the rotation rate is assumed to be much smaller than the fre- 
quency being calculated. This keeps the rotational perturbation small so that including only 
the first one or two terms in the power series expansion is satisfactory. Small is, of course 
a vague term, and it is not clear how small is small. Based on discussions at the Workshop 
on the Future of Asteroseismology held in Vienna in September 2006, estimates of the lim- 
iting rotation rate ranged from 50 to 300 km s -1 . Of course, the limiting surface equatorial 
velocity will be dependent on the mass of the star in question. 

Efforts to more accurately include rotation have been developed. These methods require 
2D calculations, so are more time consuming and c omplex. As a result, pr evious studies have 
all faced restrictions and limitations. For example, lEspinosa et al. I (120041 ) calculated the adi- 
abatic oscillations of rapidly rotating stars with uniform rotation. To succeed, they applied 
the Cowling approximation, neglect ed the Coriolis force and n eglected the Brunt- Vaisala 
frequency in the adiabatic equation. lYoshida fc Eriguchil (120011 ) have modelled quasi-radial 
modes at a range of rotation rates in rotating neutron st ars using the relativistic Cowl i ng ap - 
pro ximation. Other methods, such as t hat employed by iLignieres. Rieutord fc Reesd (120061 ) 
and lReese. Lignieres fc Rieutordl (120061 ) have fewer physical restrictions, but have so far been 
restricted to explorations of polytropic models. 

The effects discussed in this paper are only expected to matter fo r stars undergoing 
moderate to rapid rotation. A recent study of OB stars (iDaflonl 120071 ) found that 50 % 
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of OB stars have rotation velocities greater than 100 km s _1 . At least some of these stars 
are ex pected to pulsate. For example, 3 Ce phei-type pulsations have been detected in 
Spica (jSterken. Jerzykiewicz fc Manfroidl Il986l ). which is also rotating with a vsim ~ 160 
km s _1 . For the 3 Cephei stars as a categ ory, the projected rotation velocities range from 
to 300 km s _1 (IStankov fc Handler! 120051 ). The average vsim ~ 100 km s -1 , although this 
could be a selection effect, as the highest amplitude pulsators are the more slowly rotating 
stars. Another category of pulsating stars, the low amplitud e S Scuti stars (LADS) have 
been detected with vsim up to 250-300 km s -1 ( iBregerl 120071 ) . The models we consider in 
this paper are 10 M© ZAMS models with solar (Z = 0.02) metallicity. Although 3 Cephei 
stars have evolved along the main sequence, the trends produced by these models should be 
comparable to typical 3 Cephei stars. One effect which may be important is mode bumping, 
which will appear in real 3 Cephei stars, but does not appear in our unevolved models. Our 
models include uniform rotation at rates from to 0.89 fl c . Our method also allows us to 
consider differential rotation, and this will be discussed in a future paper. 



Clementl (119981 ) has developed a finite difference method for directly evaluating the 



eigenfunctions on a 2D g rid. In this paper, w e combine this method with 2D stellar mod- 
els produced by R0T0RC (jDeupred Il990l . Il995l ). The combination of these two approaches 
bypasses many of the restrictions faced by previous approaches. Our numerical methods 
and models are described in more detail in $2j We investigate the effects of rotation on the 
calculated eigenfrequencies (§3]) and eigenfunctions (Q, with the aim of establishing the 
range of validity of modes calculated with one spherical harmonic. In §5] we compare our 
results with those predicted by second order perturbation theory. 



Method 



Our st ellar models are calculated using the 2D stellar evolution code R0T0RC (IDeupree 



1990 



19951 ). allowing us to self-consistently model the surface and structure of the star 
for rotation rates from zero up to near-critical rotation. In this paper we focus on uni- 
formly rotating 10 M^Z AMS models with X =0.7, Z=0.02. We use the OPAL opacities 
Iglesias fc Rogers! (119961 ) and equation of state iRogers. Swenson fc Iglesiasl (119961 ) in these 
calculations. These models are fully 2D, with 10 angular zones from pole to equator and 349 
radial zones. We have computed a few models using 20 angular zones and find differences 
in the horizontal variation of the density to be only about 0.1%. The pulsation code uses 
Fourier transform interpolation to convert from our angular zoning to its own angular zon- 
ing, and we feel the R0T0RC angular zoning is not a major source of error in the calculations 
and use our 10 angular zone models in this work. 
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The location of the surface of the stellar model is found by assuming it lies on an 
equipotential surface. The value of the equipotential is determined by the value of the total 
potential in the angular zone which has the largest radius (for uniformly rotating models, 
this is always at the equator). The radial zone which has this value of the total potential is 
found at each angular zone and the surface boundary conditions applied there. One source 
of inaccuracy is that a radial zone is either completely interior or exterior to the surface, so 
that the surface is defined as the radial zone interface which is closest to the location of the 
equipotential. Our rotating models are made by imposing a surface equatorial velocity and 
an internal angular momentum distribution (in this case, uniform rotation) and allowing the 
surface to change as needed. This can lead to small differences between the imposed (target) 
surface equatorial velocity and the actual surface equatorial velocity, typically less than 2 
km s _1 . Throughout this paper, we refer to models by the target surface equatorial velocity. 

For our pulsa tion calculations, we use the non-radial oscillation code (NRO) developed 



by IClementl (119981 ). Instead of expressing the solution as a sum of spherical harmonics, the 
code solves the perturbation equations on a 2D spherical grid. In ROTORC, the stellar model 
is defined on a spherical polar grid, with the stellar surface location being an equipotential 
surface as discussed above. NRO tranforms this into a model defined on surfaces of constant 
density. The 2D nature of the code allows us to account for the effect of the centrifugal 
distortion, but the Coriolis force is neglected. The pressure perturbation can be expressed 
in two ways: 

oo 

SP(r,e,(f>;l,m,n) = e im ^ a™(r; n, /)P™(cos#) 

l=m 

or 



e im4, A%(r,9; n, l)r k . 



k—m 



In this code, the second form of this general equation is used. Keeping this general solution 
in mind, the linear adiabatic pulsation equations can be recast using 5 variables, related to 
the radial and angular velocity perturbation, the pressure and gravity perturbations, and the 
radial derivative of the gravitational perturbation. These variables are defined as follows: 



r k ~ 1 sin m 6 ' 
He 



v% = 



r k ^ 1 sin m ~ 1 8cos6 ' 



Us = h ■ ma' 
r k sm m u 
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Vi r k sin m Q'' 
and y 5 = <9 r y 4 

where is the radial exponent, m is the azimuthal quantum number, and k = and m = 
are special cases. If k-1 and m-1 are negative, they are replaced by 1. This form of the 
equations allows the boundary conditions to be applied while avoiding singularities. With 
these variables, the relevant linearized equations can be expressed in the general form 

drVi = f(Vh drVj^i, deVi) (4) 



The full form of the equations and their derivations can be found in IClementl f 19981 ). 



The coefficients of the finite difference expressions of the equations (as represented in 
Eqn. H]) covering the entire 2D grid can be put in a band diagonal matrix. Each element of 
this band diagonal matrix is itself a matrix, containing the coefficients at each zone in the 
2D grid. The solution of the finite difference pulsation equations proceeds in two steps, from 
the center outwards and from the surface inwards. Each integration also requires an initial 
guess of the eigenfrequency. 

The inward and outward integrations of the eigenf unctions are required to be continu- 
ous at some intermediate fitting surface. Once all of the coefficients of the equations have 
been evaluated, a subset of the matrix, including the fitting surface and the radial zones 
immediately surrounding it can be inverted to solve for the perturbations at the fitting sur- 
face and the radial zone either directly above or below the fitting surface. These values can 
then be used to step inwards and outwards through the mesh to solve for the perturbations 
throughout the rest of the grid. At some point on the surface, one of the perturbations is 
forced to be a constant (typically, 5r/r — 1) to eliminate the trivial solution of all variables 
being zero everywhere. As a result of this, there is one condition that has not been used. 
This can be used to evaluate a discriminant, which will only be satisfied (equal to zero) if 
an eigenvalue has been located. Using this method, we can step through eigenfrequency 
space, solving the matrix, evaluating the discriminant and looking for zero crossings. Once 
a crossing has been located, various convergence schemes can be used to calculate the exact 
eigenfrequency. This method can miss frequencies when two eigenfrequencies are quite close 
together, although these can usually be avoided by reducing the frequency step size. 

The code can include up to nine angular zones in the solution for the eigenfunctions, 
performing one radial integration for each angle included. At the end of the calculation, the 
solution is known at N angles, which can subsequently be decomposed into the contributions 
of individual spherical harmonics. This is done with Fourier transforms, which transform 
the N discrete points into coefficients of the appropriate cosine series. After some algebraic 
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manipulation, this series is converted into a Legendre series, which gives us the relative 
contribution of each Y™(or Legendre Polynomial for the case where m — 0). Because each 
radial integration contains angular derivatives, also evaluated using finite differencing, the 
resultant coupling among spherical harmonics arises naturally. Thus, this method allows us 
to directly model the coupling among spherical harmonics in a single pulsation mode for 
rotating stars in a natural way. 

Because I is not a legitimate quantum number for rotating models, specifying I is not 
necessary. In the pulsation code the input value of / is used to specify the parity of the mode, 
not the exact value of /. Based on the parity of I, the code includes the first k even or odd 
basis functions, where k is the input value of the number of angular zones to be included. 
We limit ourselves to small input values of / because those are expected to be the most easily 
observable. We also restrict ourselves to axisymmetric modes (m=0), although this is not 
a constraint intrinsic to the method. We have also restricted ourselves to modes with small 
radial quantum number (n). 

Because / is no longer a valid quantum number, we need a new designation for mode 
identification. We have chosen to identify the mode with a quantum number, l , which 
is the value of I of the mode in the nonrotating model to which a given mode can be 
traced back. This tracing back is based on examining both the eigenfrequency and the 
angular shape of the eigenfunction (the modes at different radial quantum number are easy 
to resolve; no mode bumping is exhibited in these ZAMS models). This is quite easy up to 
moderate rotation rates because one spherical harmonic tends to dominate. This method 
fails for rotation velocities above 420 km s -1 because no spherical harmonic dominates. For 
rotation velocities above 360 km s -1 , we find this method becomes somewhat uncertain 
and produces an irregular progression in frequency for some modes. We thus consider the 
pulsation properties for models up to 420 km s" 1 , but regard the frequencies above 360 km 
s _1 as uncertain. Although we only consider pulsation up to 420 km s -1 , our static models 
go up to near critical rotation. 



3. Accuracy of Eigenfrequencies 

As described in the above section, NRO, combined with 2D structure models from 
R0T0RC, allows us to calculate the eigenfrequencies for a rotating star without making any a 
priori assumptions about the structure of the star. The method of solution of NRO allows 
for the inclusion of multiple spherical harmonics. As a result, we can calculate eigenfunctions 
for distorted stars including the coupling between spherical harmonics. In contrast most cur- 
rent calculations and observations generally assume that pulsation frequencies and observed 
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modes can be characterized by a single spherical harmonic. It is therefore of interest to 
determine at what surface equatorial velocity modes can no longer be adequately described 
by a single spherical harmonic. 

One of the issues arising out of the following discussion is where a difference between 
two calculated modes becomes significant. Both ground-based and space-based observations 
continue to improve, as new projects are continuously launched (figuratively and literally). 
As an example, COROT is expected to measure frequencies to a precision of less than 
0.01/xHz for the lon g runs, and better than 0.065/iHz for a faint object during short runs 



(IMichel et al. 1120071 ). Based on these numbers, calculated frequencies do not need to change 
by much to be outside the observational uncertainties. However, we must ask ourselves 
whether it is reasonable to expect our models to match this accuracy. The linear adiabatic 
pulsation code uses 10~ 6 as the convergence criterion on the discriminant described in §2J 
There will be other sources of error on the final eigenfrequency, such as from the finite 
difference representation of the pulsation equations. Neglecting these other sources of error, 
NRO converges modes to an accuracy of about 10~ 6 , or about 0.001 /zHz, more than sufficient 
to match the predicted COROT accuracy. However, there are inaccuracies that result from 
the finite difference zoning in the static models. When we change the surface equatorial 
velocity from one model to the next, we change the distribution of material in the star, 
although the radial zoning (fractional surface equatorial radius) remains the same. The 
changes become larger as the rotation rate increases. This is equivalent to changing the 
radial zoning, which experience from the early calculations of linear radial pulsation indicated 
a sensitivity on about the one percent level. We have also fairly dramatically rezoned a 
couple of our models and found that the eigenfrequencies changed on about the one percent 
level, or about 8.5 /iHz for our models. The higher radial order p modes are slightly more 
affected because the outer layers of the model, where the gradients of model quantities are 
steeper, play a larger role. Clearly, our ability to measure observational frequencies to high 
precision is irrelevant until models improve enough to match them. Until then, for changes 
induced by rotational effects to be considered significant, they must be larger than our model 
uncertainties. 

Another uncertainty consideration is the angular resolution of our pulsation calculations. 
As described above, the number of spherical harmonics used in NRO determines the number 
of radial integrations performed. There are several ways we can assess the effects of this 
changing angular resolution. Firstly, we would expect the slowly rotating modes to be 
relatively unaffected by angular resolution. This is indeed what we find. In the case of slow 
rotation, the coefficients for the higher order spherical harmonics are small, typically not 
more than a percent up to 120 km s _1 . Over these same rotation ranges, we also expect the 
frequency to be relatively unaffected by angular resolution, and this is indeed what we find. 
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The frequencies shown in Fig. [T] differ by less than a quarter of a percent over this rotation 
range. 

In the majority of our plots, we show our results as a function of surface equatorial 
velocity, as this is the unit most easily compared to observations. However, for comparison 
with other models, it is more useful to show results as a function of angular rotation rate 
expressed as a fraction of critical rotation (Q/Q c ). Critical rotation was calculated using a 
model rotating at 575 km s _1 , with an equatorial radius of 5.792 R . This model is quite 
close to critical rotation. We have summarized the conversion between these two frames of 
reference, as well as some other parameters of our models in Table [TJ 

3.1. Frequency Changes 

The simplest way to determine where the assumption that a single Y™can be used is to 
compare the frequencies as calculated with different numbers of spherical harmonics. This 
is illustrated in Fig. [TJ which shows the normalized frequencies for the Iq = and Iq = 1 
fundamental modes, as calculated using 1, 2, 3 and 6 spherical harmonics. At some cut off 
surface equatorial velocity, the eigenf unctions calculated with only a few spherical harmonics 
begin to deviate significantly from those calculated using 6 spherical harmonics. For the Iq 
= mode, the frequencies calculated with 1 spherical harmonic are in reasonably good 
agreement to quite high velocities, remaining within approximately 0.5 % of the frequencies 
calculated with more spherical harmonics. The Iq = 1 mode as calculated with 1 spherical 
harmonic rapidly diverges from the frequencies as calculated with multiple basis functions. 
In this case, the single spherical harmonic frequency reaches a difference of 1% at a surface 
equatorial velocity of 210 km s _1 (0.51 Q c ). 

Similar results are found for higher order modes. These results are summarized in Table 
[2j To determine the location of the cut off surface equatorial velocity, as described above, 
we take a difference of 1% to be significant, as discussed in SJ21 

Although the periods are expected to change depending on the details of the model, 
period differences and ratios are expected to be much more stable. Hence, in the next two 
sections we will consider the large separation and period ratios of our frequency calculations. 
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Table 1. Summary of model parameters 
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Fig. 1. — The frequency changes as a function of surface equatorial velocity for the funda- 
mental mode for l = (top) and Iq = 1 (bottom). Frequencies shown are calculated with 
(O) - 1 Y™, o - 2 Y™, (□) - 3 Y™ and (A) - 6 Y™. 
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3.2. Large Separations 

We have studied the large separation between the n — 0, 1 and 2 modes for l = 0-3. 
We have calculated the large separations in the usual way 

Av = vi tn+1 - Vt >n . (5) 

Before comparing these for the effects of the number of spherical harmonics included, we 
need to account for rotation, which can change the large separation by changing the model 
structure. First, we normalize these large separations with respect to the non-rotating model 

Dv = Au(v = 0) - Au(v). (6) 

We can then use these normalized large separations to look for the effects of the number of 
spherical harmonics included in the calculation (JV) 

Vv = Du N — Dv N=6 . (7) 

For this calculation, we have assumed that the frequencies calculated with 6 Y™s are closest 
to the true pulsation frequencies, so the smaller the differences between this and other 
calculations, the more accurate the smaller number of spherical harmonics. This is illustrated 
in Fig. [21 which shows the results of Eqn. [7] as a function of surface equatorial velocity for 
the separation between the Iq — fundamental and first harmonic. 

The uncertainty in the theoretical calculations of large separation is inversely propor- 
tional to the uncertainty in the radius of the stellar model in question. Taking the uncertainty 
in radius to be the size of one radial zone, for our models, this is approximately 0.04yuHz. 
Observationally, large separations are well determined for solar type stars, with uncertainties 
typically less than 1/iHz. As a conservative estimate, we have chosen 1/zHz as our significance 
criterion, as shown by the dashed lines in Fig. [2J It should be noted that once the large 
separations with 1 and 2 spherical harmonics begin to diverge, they do so quite rapidly, so 
unless the cut off criterion is appreciably smaller ( < 0.5/iHz), the cut off surface equatorial 
velocity is not an extremely sensitive function of the cut off criterion. The limiting rotational 
velocities estimated using the large separations are summarized in column 4 of Table [2J 



4. Accuracy of Eigenfunctions 

So far, the limiting rotation rates entered in Table [2] have been for the lo = and 1 modes 
only. This is a result of the way spherical harmonics are included in NRO. To calculate the lo 
= 2 mode, for example, the code will select even Y™s starting with l = 0, so at least 2 Y£"s 
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Equatorial Velocity (km/s) 



Fig. 2. — The relative large separation (Eqn. [7j) as a function of surface equatorial velocity 
between the Iq = fundamental and first harmonic. Symbols are as follows: (O) - 1 spherical 
harmonic, (A) - 2 spherical harmonics (o) - 3 spherical harmonics, all relative to 6 spherical 
harmonics. Dashed lines indicate the significance criterion adopted in this work. 
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are required. This is true for any mode with Iq > 2. As a result, we cannot directly compare 
eigenfrequencies calculated with several spherical harmonics to those calculated with a single 
spherical harmonic. 

We can still compare the eigenf unctions, and in this section this is what we will do. One 
of the advantages of including several spherical harmonics is the ability to study the effect 
of rotation not only on the eigenfrequencies, but also on the shapes of the eigenf unctions. 
For a non-rotating object, regardless of how many spherical harmonics are included, the 
eigenfunction remains a pure Y™,as it should. As the rotation rate increases, neighbouring 
spherical harmonics begin to contribute progressively more to the shape of the eigenfunction. 
These effects could be quite important for mode identification, and need to be considered in 
rapidly rotating stars. One technique for mode identification uses the pulsation amplitudes 
in different colors as determined by single spherical harmonics. With rotation significantly 
altering the modes by coupling spherical harmonics, it could alter these color amplitudes 
and change the mode identification. We find that the effects of the coupling can become 
significant, even at very moderate rotation rates. 

We have used a combination of the value of the eigenfrequency and the angular variation 
of the eigenfunction at the surface to identify the modes as we progressed from one rotating 
model to the next. Of course, with the finite difference approach the angular variation of the 
eigenfunction can vary with depth. Fig. [3] presents this variation for several rotation rates for 
the radial fundamental mode. Each plot contains the variation at several different depths. 
As expected, the variation with depth is small for slowly rotating models, and grows as the 
rotation rate increases. Despite this growth in variation, the profile remains recognizably the 
same until the most rapid rotation rate presented. This occurs at a rotation rate at which 
we are already beginning to have trouble tracing the modes from one rotation rate to the 
next as we have previously mentioned. 

Fig. H] shows the angular variation at the surface in the radial component of the Iq = 
fundamental mode at 90 and 270 km s _1 . At 90 km s -1 , the distorting effects of rotation 
are negligible, although the differences are visible. In contrast, by 270 km s _1 the differences 
between the numbers of spherical harmonics are quite significant, and 1 spherical harmonic 
is clearly not sufficient to model the horizontal shape of the mode. In comparison, the 
eigenfrequencies were considered to be accurate using one spherical harmonic up to rotation 
rates of 300 km s _1 . This highlights the truism that even marginal eigenf unctions can give 
reasonable eigenfrequencies. By 270 km s _1 , the mode no longer looks like an I = mode, nor 
even an Z = 2, but is beginning to distinctly show the characteristics of the I = 4 contribution. 
These two velocities were chosen based on the relative contribution of each Yj 11 , shown for 
the radial fundamental mode in Fig. [51 At 90 km s -1 , with all three sets of basis functions, 
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Fig. 3. — Variation in the radial eigenfunction for the l Q mode as a function of colatitude 
at various depths (fractional surface equatorial radii of approximately 0.1, 0.2, 0.3, 0.4, 0.5, 
0.6, 0.7, 0.8, 0.9 and 1.0) for models rotating at 50, 150, 240, 300, 360, and 420 km s _1 . 
The convective boundary is located between 0.2 and 0.3 K eq . The variation at each depth is 
normalized to be unity at the pole for purposes of comparison. The variation is smallest at 
the center of the star, and increases towards the surface. On the plot of the 420 km s" 1 , the 
layer closest to the center is indicated with a dashed line, and the layer closest to the surface 
is indicated by a dot-dashed line. In most cases, 420 km s _1 is the most rapidly rotating 
model considered, as mode identification becomes difficult. 
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Fig. 4. — Angular variation in the radial eigenfunction for the radial fundamental mode 
of a model rotating at 90 (top) and 270 km s -1 (bottom). On both plots, the shape of the 
eigenfunction is shown as calculated using 1 (dotted), 2 (solid), 3 (dashed) and 6 (dot-dashed) 
spherical harmonics. 
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the I = component contributes nearly 100%, while at 270 km s l , the contribution of the 
same component drops below 50% when 6 spherical harmonics are considered. 

From Fig. we can see that with two and three spherical harmonics, all of the spherical 
harmonics contribute a relatively significant amount by the time the model is rotating at 
intermediate speeds. In contrast, with six spherical harmonics, the contribution from the 
highest order spherical harmonics (I = 10) remains small out to at least 300 km s _1 . Although 
the contribution starts to become significant at very high rotation rates (v > 350km s _1 ), it 
still remains a factor of 2-3 lower than the main contributors. From this, we have taken the 
shape of the eigenfunction with six spherical harmonics as being the most correct and have 
used it as a basis of comparison. 

Based on the results shown in Fig. HJ we know that one spherical harmonic ceases to be 
sufficient somewhere between 90 and 270 km s -1 . From Fig. [5J we can see that the relative 
contribution of the F m drops below 90% at a surface equatorial velocity between 150 and 
180 km s _1 . The angular variation of the eigenf unctions for these two velocities are shown 
in Fig. O It is at this point that we would say multiple spherical harmonics are required to 
accurately reproduce the shape of the mode (cf. Fig. E]). 

We have developed a quantitative measure of how the shapes of the eigenfunction differ 
from that calculated using six spherical harmonics. This estimate is calculated by taking the 
absolute value of the difference between the 6 basis function eigenfunction (standard) and 
one of the other eigenfunctions (comparison) at 9 points. These points are equally spaced 
across the surface of the model, with 9 = lOi. The point at 9 = is excluded, as all the 
eigenfunctions are normalized to one at this point. These differences are then squared and 
summed. The square root of the sum is normalized by the number of points to give a measure 
of how different the two curves are: 



This difference as a function of surface equatorial velocity is shown in Fig. [7J The differences 
between the eigenfunctions calculated with 1, 2 and 3 spherical harmonics relative to 6 
spherical harmonics rises sharply starting at a surface equatorial velocity of 180 km s _1 . 
Based on this rise and the eigenfunctions shown in Fig. [6j we estimate that when the mean 
difference rises above 0.06, more spherical harmonics are needed to accurately reproduce the 
shape of the mode. 

For the other modes, the results are qualitatively similar, although the extent of the 
differences varies. The results for all four l values considered in this paper are summarized 
in Table [2j Overall, one spherical harmonic remains a good approximation out to at least 
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Equatorial Velocity (km/s) 

Fig. 5. — The relative contribution to the F mode of each spherical harmonic for 2 (top) 
3 (middle) and 6 (bottom) spherical harmonics. In the top plot, after v ~ 150 km s _1 , 
the contribution from Iq = drops below ~ 90% and we say that you need more spherical 
harmonics to be able to model the mode. Symbols are defined as follows: O - I — 0, □ - 
I = 2, x - I = 4, o - I = 6, + - I = 8, A - I = 10. 
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Fig. 6. — As for Fig. HI but for the velocities on either side of the cut off surface equatorial 
velocity. At the lower velocity (150 km s -1 , top), the shape can be calculated reasonably 
well using one Y[™, but at the higher velocity (180 km s _1 , bottom), 2 or more are needed 
to accurately reproduce the horizontal variation in the eigenfunction. Symbols are the same 
as in Fig. HI 
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Equatorial Velocity (km/s) 

Fig. 7. — The mean difference between the shape of the radial fundamental eigenfunction 
with 6 spherical harmonics and a pure P mode (O), 2 spherical harmonics (□) and 3 spherical 
harmonics (A). Although there is some variation, all three curves show a sharp rise beyond 
200 km s _1 . See text for the definition of the mean difference. 
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90 km s _1 (0.23 Q c ). For some modes, such as the radial fundamental, this approximation 
remains valid to much higher rotation rates (270 km s _1 , 0.64 Q c ). As both the angular and 
radial order of the mode increase, the limiting surface equatorial velocity decreases. In most 
cases, we find that the differences among calculations with different numbers of spherical 
harmonics grow quickly as a function of surface equatorial velocity once the differences 
become sizeable. We can conclude that our results are not particularly sensitive to the exact 
value of the cutoff criterion we have chosen, as long as it is not significantly lower than what 
we have used. We also find that comparing frequencies or frequency differences produce 
approximately the same results. Based on our results for a 10 M Q model, a single spherical 
harmonic is never a good approximation for rotation rates above 0.64 fl c , appears to always 
be a good approximation for rotation rates below 0.23 f2 c , and must be used with caution 
for rotation rates between these two values. Although there may be some mass effects, we 
do not expect these results to change significantly for masses close to 10 M . 



5. Comparison with Perturbation theory 

Second order perturbation theory is routinely used to compute linear pulsation modes for 
rotating stars in which the centrifugal forces are expected to affect the pulsation frequencies. 
It has been difficult to comment on when second order perturbation theory can be expected 
to fail because there have been few calculations of eigenfrequencies using other methods. 
Our approach will allow placing some limits on the range of applicability of second order 
perturbation theory, but again these limits will be a product of the accuracy obtainable or 
required. 

Second order perturbation theory shows that, for axisymmetric modes as we consider 
here, the change in eigenfrequency is a linear function of the square of the rotation rate (e.g., 
Saio 1981). We shall compare our results with this linear relation in two separate ways, both 
of which determine the failure of perturbation theory by a deviation from this linear relation. 
Of course, the result will depend on the quantitative value as to when the deviation becomes 
significant, a point we will discuss at the end of this section. We shall use the results we 
feel most accurately reflect the true values of the pulsation frequencies, the results with six 
angular zones in the 2D pulsation grid for our comparison of eigenfrequencies. 

The first method starts with the first four models in the rotation sequence (surface 
equatorial rotation velocities from to 90 km s~ 1 . We calculate the best fit to the linear 
relationship as given by perturbation theory, and the standard deviation. We repeat this 
exercise, each time adding one more model to the analysis, until all rotation velocities are 
included. As long as the linear relation is satisfied, we expect the standard deviation to be 
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Table 2. Summary of velocities at which 1 Y^fails to accurately reproduce the mode. 



k 


n 


Frequency 21 




Eigenf unction 








>360 




165 





1 


240 


160 


60 





2 


180 


24 


25 


1 





210 




110 


1 


1 


210 


140 


105 


1 


2 


180 


30 


85 


2 









75 


2 


1 






60 


2 


2 






45 


3 









70 


3 


1 






85 


3 


2 









a Limiting surface equatorial velocity based on frequency differences larger than 1% 

b Limiting surface equatorial velocity based on difference in the large separation greater 
than 1/iHz 

c Limiting surface equatorial velocity based on eigenfunctions with mean differences larger 
than 0.06 
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approximately constant as we add results for more rapidly rotating models. At some point, 
as the rotation becomes more rapid, the standard deviation will become larger and at some 
threshold value will be declared no longer to be an adequate representation of a straight 
line. Thus second order perturbation theory would no longer be considered reliable. We 
plot this standard deviation as a function of the rotation rate of the most rapidly rotating 
member of each sample in Fig. [HJ We somewhat arbitrarily set our threshold at 4 x 10~ 6 as 
being a value above the flat region for all modes. The values for the limits of applicability of 
perturbation theory computed by this method are listed in the column entitled 'Linear Fit' 
of Table [3j We have also examined the slope of each linear fit, and as expected, find that 
the slope changes gradually where the linear fit is good, and more rapidly as more points 
are added. 

One difficulty with the above approach is that the coefficients of the linear fit change as 
more rapidly rotating models are added. A more constraining determination of the threshold 
of perturbation theory might be obtained by using the first few members of the sequence to 
determine the coefficient of the linear fit. The assumption is that the slope that perturbation 
theory would predict is correctly computed using the first few slowly rotating members of the 
sequence. We use the first five members in our rotation sequence to calculate this coefficient. 
We then use this coefficient to determine perturbation theory frequencies at each of our 
surface equatorial velocities. As before, we take the differences between the two methods 
as significant when they are larger than 1%. The results for this method are listed in the 
column entitled 'Coefficient Fit' of TableEJ We compare our pulsation frequencies with those 
predicted assuming the coefficient computed for the first four members of the sequence is 
valid at all rotational velocities in Fig. [9j 

We find the trends for both methods of evaluating the threshold are similar for the two 
methods, but that the thresholds computed for the coefficient fit are more constrained. This 
is to be expected because forcing a linear fit to have a certain slope is more confining that 
merely forcing a fit to be linear. It is interesting that the threshold for perturbation theory 
occurs at generally higher rotation speeds than the threshold for the validity of a single 
spherical harmonic. The extrapolation of the linear fit to higher rotation velocities is flatter 
than our calculation with six angular zones and much flatter than our calculation with only 
one angular zone. 

Our results indicate that per turbation theory is satisfacto r y to a ppreciably larger rota- 



tion velocities than the results of iReese. Lignieres fc Rieutordl (120061 ) . who found that third 



order perturbation theory failed for rotation rates above about 0.2 Q c . Much of this differ- 
ence arises from the much tighter constraint they placed on what difference in eigenvalues 
is significant. They are able to do this because they perform their comparisons using poly- 
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tropes, which can be numerically integrated very accurately, whereas we use finite difference 
techniques to generate our more realistic stellar models. A subsidiary consideration is that 
they can control both the total mass and radius, and thus can arbitrarily scale from one 
model to the next, whereas our models include the conservation of energy, which removes 
the radius as an arbitrary parameter. Also, the surface locations at each angle of our rotat- 
ing models are quantized; the surface is regarded to include the full radial zone instead of 
fractions of zones. Our errors are in line with variations in eigenvalues computed for radial 
modes at a similar stage of development (e.g., Castor 1971). We believe these errors are 
reasonable at the present time because the deduced properties of the stars observed will 
be inaccurate both from the conversion from observed parameters to theoretical parameters 
and from the uncertainties in the effects of inclination on the relation between the observed 
and intrinsic properties. The model and parameter inaccuracies will be far greater than the 
error in the observed frequencies. Physical uncertainties, particularly in the internal angular 
momentum distribution, are expected to be greater than or equal to the uncertainties in an 
individual model, particularly for the more rapidly rotating stars in which we are interested 
(v i 200km s _1 ). We believe that being able to compute the evolution of the rotation law 
as the star ages may, at this stage, play a more important role than increasing the accuracy 
of the calculations. Of course, we recognize that improvements in accuracy on all fronts are 
valuable. 



6. Conclusions 

In this paper, we have attempted to test the validity of two independent assumptions 
commonly made in calculating stellar oscillation frequencies. These are firstly, that the 
non-radial modes can be modelled using a single Y™, and secondly, that the modes can be 
calculated using second order perturbation theory out to some limiting (highly uncertain) 
rotation rate. 

We find that when a single spherical harmonic becomes inaccurate is mode dependent, 
with it failing at lower rotation velocities for higher order modes. The answer is also different 
depending on what property one examines. A single spherical harmonic is sufficient to 
reproduce frequencies to within 1% for rotation velocities up to at least 180 km s _1 (0.44Q c ), 
and for some low order modes, may even be valid up to 390km s _1 (0.85Q c ). In contrast, the 
angular shapes of the eigenfunctions are extremely sensitive to rotation, and the assumption 
fails at a maximum surface equatorial velocity of 165 km s _1 . In most cases, the assumption 
fails at much lower rotation velocities, typically around 50- 75km s _1 . Period differences (large 
separations) are expected to be of most interest, and these are also found to be sensitive to the 
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Fig. 8. — Standard deviation from a straight line as more points are included for the Iq = 
and 1 modes (top) and Iq — 2 and 3 modes (bottom). We take the cut off standard deviation 
to be 4xl0 -6 . Symbols are as follows: O - fundamental, □ - first harmonic, A - second 
harmonic. Solid lines represent the even modes (0, 2) and dashed lines represent the odd 
modes (1, 3). 
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Fig. 9. — Normalized frequencies as calculated with NRO (o) and using an estimate of the 
perturbation theory results (x) for the l — 2 f mode. 



-28- 



Table 3. Summary of velocities at which perturbation theory fails to accurately reproduce 

the mode. 
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order of the mode. A single spherical harmonic can accurately predict the difference between 
the fundamental and first harmonic of the lo = and 1 mode up to velocities of around 150km 
s _1 (0.37 Q c ). The higher order modes are very sensitive to rotation, and the assumption fails 
at velocities of around 25-30 km s _1 (0.08 fl c ). One interesting consequence of the limitations 
of a single spherical harmonic is the impact it may have on mode identification, which is 
most often based on comparing the variation in pulsation amplitude with color with models 
computed assuming a single spherical harmonic (e.g., Heynderickx, Waelkens & Smeyers 
1994). 

We have compared our eigenfrequencies with the relation between eigenfrequency and 
rotation rate predicted by second order perturbation theory. The relationship is followed 
reasonably well for models rotating up to surface rotational velocities of about 400 km s _1 for 
very low order modes. The relation fails at lower rotational velocities (approximately 200 km 
s _1 or f2/f2 c = 0.58) for modes with two or three radial nodes. These values are dependent on 
the difference between the two sets of frequencies tolerated. In these calculations, the limits 
are determined by the properties of the rotating stellar models rather than the calculations 
of the eigenfunctions. 

The authors wish to thank M.J. Clement for the use of NRO and for his assistance in 
using and understanding the code. This research was supported by a Natural Science and 
Engineering Council of Canada (NSERC) Discovery grant and a NSERC graduate scholar- 
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Innovation and the Nova Scotia Innovation Research Trust. 
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